﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Nature;
using Nature.Numerics;
using System.Diagnostics;
using Nature.Numerics.Basic;

namespace Nature.Toolkit.GEqu
{
    public class TVIdealGasChemicalEquilibriumState : IChemicalEquilibriumState, IFunction
    {
        #region Private Fields
        IChemicalEquilibriumState m_gibbs;
        bool m_equilibrated = false;
        double[] m_x0;
        /// <summary>
        /// mixture mass density
        /// </summary>
        double m_ro;
        double? m_t0;
        double? m_p, m_p0;
        bool m_compositionInitialized = false;
        Newton<TVIdealGasChemicalEquilibriumState> m_newton;
        #endregion

        #region Private Members
        void Equilibrate()
        {
            if (m_equilibrated) { return; }
            if (false == m_t0.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullTemperatureException(); }
            if (false == m_p0.HasValue) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullPressureException(); }
            if (false == m_compositionInitialized) { throw ChemicalEquilibriumStateNotInitializedException.CreateNullCompositionException(); }
            double c0 = (m_p0.Value / m_t0.Value / Constants.Rgas); // initial mixture concentration [mol/m**3]
            double w0 = this.Model.CalcMixtureMolarMassFromX(m_x0); // initial mixture molar mass [kg/mol]
            m_ro = c0 * w0; // fixed mixture mass density
            m_p = m_newton.Solve(m_p ?? m_p0.Value);
            m_equilibrated = true;
        }
        #endregion

        public TVIdealGasChemicalEquilibriumState(ModelObject model)
        {
            this.Model = model;
            m_gibbs = new TPIdealGasChemicalEquilibriumState(model);
            m_newton = new Newton<TVIdealGasChemicalEquilibriumState>(this);
            m_x0 = new double[model.Nsp];
        }

        public ModelObject Model { get; private set; }

        public int Nsp { get { return this.Model.Nsp; } }

        #region IChEqCalculator Members

        double IChemicalEquilibriumState.T
        {
            get
            {
                if (m_t0.HasValue) { return m_t0.Value; }
                throw ChemicalEquilibriumStateNotInitializedException.CreateNullTemperatureException();
            }
            set
            {
                if (m_t0 != value)
                {
                    m_equilibrated = false;
                    m_t0 = value;
                    m_gibbs.T = value;
                }
            }
        }

        double IChemicalEquilibriumState.P
        {
            get
            {
                Equilibrate();
                return m_p.Value;
            }
            set
            {
                if (m_p0 != value)
                {
                    m_equilibrated = false;
                    m_p0 = value;
                }
            }
        }

        ReadOnlyArray<double> IChemicalEquilibriumState.X
        {
            get
            {
                Equilibrate();
                return m_gibbs.X;
            }
            set
            {
                m_compositionInitialized = true;
                bool changing = false;
                for (int isp = 0; isp < Nsp; ++isp)
                {
                    changing |= (m_x0[isp] != value[isp]);
                    m_x0[isp] = value[isp];
                }
                if (changing)
                {
                    m_equilibrated = false;
                    m_gibbs.X = value;
                }
            }
        }

        void IChemicalEquilibriumState.DropState()
        {
            m_p = null;
            m_t0 = null;
            m_p0 = null;
            m_compositionInitialized = false;
            m_equilibrated = false;
            m_gibbs.DropState();
        }

        #endregion

        #region IFunction Members

        double IFunction.Compute(double pressure)
        {
            m_gibbs.P = pressure;
            double w = this.Model.CalcMixtureMolarMassFromX(m_gibbs.X); // mixture molar mass [kg/mol]
            double c = (m_gibbs.P / m_gibbs.T / Constants.Rgas); // mixture molar concentration [mol/m**3]
            double diff = m_ro - c * w;
            return diff;
        }

        #endregion
    }
}